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I. INTRODUCTION 



In this paper we present our final results for the renormalization and improvement con- 
stants for quark bilinear operators using Wilson's gauge action and the 0{a) improved Dirac 
action first proposed by Sheikholeslami and Wohlert ■1]. The calculations have been done 
at three values of the gauge coupling, /3 = 6.0, 6.2, and 6.4 in the quenched approximation.^ 
Our results represent a realization of Symanzik's improvement program for systematically 
reducing discretization errors in lattice simulations |3i,i^]. Results for the improvement of 
the Dirac action have been obtained previously by the ALPHA collaboration and we have 
used these in our calculation. This paper deals with the improvement of external bilinear 
operators, O with O being one of the five Lorentz structures A, V, P, S, T. 

The mixing with extra operators, both for on-shell and off-shell improvement of the 
operators, and the introduction of mass dependence in the renormalizaton constants has 
been discussed in detail in Section II of Ref. 0. To summarize that discussion, and to 
remind the reader of the notation, the fully improved and renormalized bilinear operators 
at 0{a) are 

(Af)^ = Z° (1 +bAamj) {A^ + acAd^P - a^c'^(V^^*^^75>V^^''^ - 

(TP),, = Z^il + bTam,,) {T,, + acrid^V, - d,V,) - a^^^i^^'ka^^W^^^^ - ^^'^ma,,^^^^)) 
(PP) = zUl + bpamij){P - ayp{lp^'\,W^P^^^ -i^^'^W^^ 
(5f) ^ Z'si^ + bsam,,){S-a^c's{^^'^W^^^^-i^^^^^ 

Here {ij) (with i ^ j) specifies the flavor. The Z^ are renormalization constants in the 
chiral limit and rhij is the quark mass defined in Eq. ((7j) using the axial Ward identity 

(AWI). Wipj = + + 0{a?) is defined to be the full 0(a) improved Dirac operator 

for quark flavor j (See Appendix in Ref. 0). This ensures that the equation of motion 
operators give rise only to contact terms, and thus cannot change the overall normalization 
Zo- The normalization is chosen such that, at tree level, = 1 for all Dirac structures. 

We determine the improvement and renormalization constants using Ward identities. 
When implementing these, we have a number of choices. Two are of particular importance. 
First, we need to pick a discretization of the total derivatives appearing in the improvement 
terms proportional to CAy,T- Note that, because the derivatives are external to the operators, 
rather than internal, this choice should not impact the result for the Z's or c's, aside from 
corrections of 0{a?) which are not controlled. In fact, we will find that such higher order 
corrections are largely kinematical, and can be removed by the chiral extrapolations. Second, 
we need to choose the external states. As far as we know, there are no standard choices, and 
so we take either the state giving the best signal, or an average if there are several giving 
similar accuracy. We then use the difference of the results with those from other states as 
part of the estimate of the uncertainty. Although this is somewhat ad hoc, it is a well-defined 



Preliminary results were presented in |2| and are updated here. 
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procedure as long as we make consistent choices for all lattice spacings. We stress that the 
coefficients bo differ from the bo used by earlier authors. These are related as 

boaniij = boafhij + 0(a) (1) 

where rriij = {rrii + mj)/2 is the average bare quark mass defined as arrii = l/2Ki — 1/2kc, 
K being the hopping parameter in the Sheikholeslami-Wohlert action and Kc its value in the 
chiral limit. At the level of 0{a) improvement, one has 

bo = {ZlZllZDbo . (2) 

The analogous relation between m and m is given in Eq. ((201) • 

In this paper we present results for those overall normalization constants, Zq, that are 
scale independent and the improvement constants bru co, and co- A detailed discussion 
of the methods has already been presented in Refs. y and 0, and we do not repeat them 
here. The extension of the method to full QCD has been presented in 0. Instead we 
concentrate on presenting the final results and new aspects of the analyses. In particular, 
using three lattice spacings we are able to significantly improve our understanding of residual 
discretization and perturbative errors by comparing our results with those obtained by the 
ALPHA collaboration using a non-perturbative method based on the Schrodinger functional 
and with the predictions of perturbation theory at one-loop order. 

The remainder of this paper is organized as follows. In the next section we describe the 
essential features of our simulations and the types of propagator we use. Section ITTTl gives 
an overview of the methods we use to implement Ward identities and a summary of the 
results. We then run through the results from the different Ward identities that are needed 
to calculate ca (Sees. II VI and IVI for zero and non-zero spatial momenta, respectively), Zy 
and by (Sec^EH), cy and 6^ -by (SecEIIl), fSec. IVTTT|l . Z^/ZO and 65 -bp (SecHXj, 
bp — bA and bs (SecEj), ct fSec.lXI|). and the coefficients of the equation-of-motion operators 
(Sec. IXII|) . We compare our results with those of others in Sec. IXIIII and with one-loop 
perturbation theory in Sec. IXIVI We close with brief conclusions in Sec. IXVI 



II. DETAILS OF SIMULATIONS 

The parameters used in the simulations at the three values of (3 are given in Table B The 
table also gives the labels used to refer to the different simulations. For the lattice scale 
a we have taken the value determined in Ref. using tq as it does not rely on the choice 
of the fermion action for a given j3. The values of the hopping parameter k, along with 
the corresponding results for the quark mass am, determined using the Axial Ward Identity 
( AWI) , and aM^ are given in Table |n] . 

Four major changes have been made in the analysis compared to our previous work j^. 
First, the addition of the data set at /3 = 6.4 to those aX (3 = 6.0 and 6.2 (the latter two being 
unchanged from Ref. |5j) allows the identification of higher order contributions in the chiral 
extrapolations. As a result we now use quadratic or linear fits in the chiral extrapolations 
for all three (3 values as opposed to the linear or constant fits used in j^. Second, the 
improvement in the signal with increasing j3 allows us to better determine which values of k 
to keep in the fits. We are able to use all seven values, ni — kj, ai (3 = 6.2 and 6.4 whereas 
(3 = 6.0 data at k = are too noisy (no clear plateaus in the ratios of correlators), and in 
some cases even the data at k = are too noisy to include in the fits. 
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Label 


/3 


csw 


(GeV) 


Volume 


L (fm) 


Confs. 




60NPf 
60NPb 


6.0 


1.769 


2.12 


16^ X 48 


1.5 


125 
112 


4-18 
27-44 


62NP 


6.2 


1.614 


2.91 


24^ X 64 


1.65 


70 
70 


6-25 
39 - 58 


64NP 


6.4 


1.526 


3.85 


24^ X 64 


1.25 


60 


8-56 



TABLE L Simulation parameters and statistics, denotes the time interval over which the chiral 
rotation is performed in the AWL The initial Wuppertal source is placed at i = 0. 

The third improvement is with respect to the discretization of the derivatives in the 
operators. As in Refs. |^ 0, we use two discretization schemes in order to estimate the size 
of O(a^) uncertainties. Most of our central values come from the "two-point scheme" (which 
is changed from Refs. H, 0). This uses two-point discretization^ throughout the calculation, 
i.e. both in the axial rotation of the action, 55*, and in the operators. It improves upon 
the scheme with the same name that we used in refs. IH S in which we only used two-point 
discretization in the calculation of ca and 5S, but all other operators were discretized using 
three-point discretization. 

We estimate discretization errors using a hybrid scheme in which we use three-point 
discretization^ in all the operators but retain the two-point discretization in 5S (using the 
corresponding two-point values for ca and m). We refer to this as the "three-point scheme". 
We did not use three-point derivatives in the discretization of 5S in the present calculation 
for reasons of computational cost. We stress, however, that both schemes have errors starting 
at 0{a?). By comparing them we obtain information about about the size of these errors. 
Further details on the two schemes are explained later. 

Lastly, we have also added the calculation of ca using a "four-point" discretization of 
derivatives^ which is improved to 0{a^) at the classical level. This allows us to further 
study discretization errors. 

The fourth improvement is in the definition of the central value x obtained from the 
jackknife fits. We now include an 0(1/A^) correction in the single elimination jackknife 
procedure and define 

X = x + N{xq — x) (3) 

where x = J^N^jk/^ is the uncorrected (and previously used) estimate, is the sample 
size, and Xq is the result of the fit to the full data sample. 

The reanalysis changes many of the results presented in The most significant changes 
(with final results changing by more than la) arise from the order and range of the fit used 



2 f{x + 0.5a)^[f{x + a) + f{x)]/2, 
dxf{x + 0.5a) — > [f{x + a) — f{x)]/a and 

dlf{x + 0.5a) ^ [/(.T + 2a) - f{x + a) ~ f{x) + f{x - a)]/{2a^). 

3 dxfix) ^ [fix + a) - fix - a)]/(2a), and 
d'Jix) ^ [fix + a)- 2fix) + fix - a)]/a\ 

4 fix + 0.5a) ^ (9[/(.T + a) + fix)] - [fix + 2a) + fix - a)])/16, 

d^fix + 0.5a) ^ (9([/(a: + a) - fix)] - [fix + 2a) - fix - a)]/27)/8a and 
dlfix + 0.5a) ^ [fix + 2a) - fix + a) - fix) + fix - a)]/i2a^). 
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(for example, changing from linear to quadratic extrapolation). The other changes in the 
analysis lead to smaller changes in the final results. We comment below on the changes at 
appropriate places. Because of these changes we present here estimates from all four sets 
of simulations listed in Table lU and these revised estimates supersede previously published 
numbers. 

To highlight the improvement in the signal in various ratios of correlation functions with 
P, we include in Figures UHl CH IHIITl 1201 121 and El previous data from 60NP and 62NP 

sets for comparison. Whereas the signal is marginal at f3 = 6.0, it improves rapidly, and 
by /? = 6.4 reliable estimates for all constants can be obtained with O(IOO) independent 
configurations. 

For each set of simulation parameters the quark propagators are calculated using Wup- 
pertal smearing The hopping parameter in the 3-dimensional Klein-Gordon equation 
used to generate the gauge-invariant smearing is set to 0.181, which gives mean squared 
smearing radii of (r/a)^ ^ 2.9, 3.9, and 5.4 for (3 = 6.0, 6.2, and 6.4 respectively. 

In Table Ewe also show the time extent of the region of chiral rotation in the three-point 
axial Ward identities. The dependence of our results on this region was investigated at 
/3 = 6.0, as shown by the two different time intervals listed under 60NPf and 60NPb. We 
observed no significant difference in the two results, so for our final results we average the 
two values weighted by their errors. In the 62NP calculation, we used two separate rotation 
regions with equal time extent and placed symmetrically about the source. This allowed us 
to average the correlation functions to improve the statistical sample. In the 64NP data 
set we were able to further improve the efficiency of the method by enlarging the region of 
insertion to include the whole lattice except for a few time slices placed symmetrically on 
either side of the source for the original propagator at t = 0. This construct allows us to 
average the signal from forward and backward propagation with a single insertion region, 
(time slices 8 — 56), and reduces the computational time significantly because only five 
inversions are required instead of the eight needed in the 60NP and 62NP studies (where 
forward and backward propagating correlators were calculated separately). 

The five kinds of propagators we use in our calculation at /5 = 6.4 are as follows. The 
initial quark propagator is calculated with a Wuppertal source on time-slice t = for all 
the lattices. To make explicit the construction of sources for propagators with insertions 
we label the two ends of the time integration region by (tj,tj), which for 64NP data are 
tj = 8 and tf = 56 as listed in Table EJ We define the insertion operator 6S using the 
two-point discretization of the derivatives, whereby the discretizations for the three terms 
in SS = 2mP — d^A^ — acAdlP are 

J d^xd^Ai J d^x[A4{tf,x) — Ai{ti,x)] , 

Jd^xdlP^ J d^x[Pitf + l)-P{tf-l)-P{U + l) + P{ti-l)]/2, (4) 

j d^xP ^ J d'x[P{U)/2 + P{U + 1) + . . . + P{tf - 1) + P(t/)/2] . 

Starting with the original Wuppertal source propagator we construct the three quantities 
defined in Eq. [Hand use these as sources to create the propagators with insertions. The final, 
fifth, propagator is calculated by inserting 75 at zero 3-momentum on time slice t = 23, 20 
and 24 respectively for the three (3 values. This is needed to study the vector Ward identity 
used to extract Zy- 
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The quark and antiquark in the operators in SS, which have flavors we call "1" and "2" 
respectively, are always taken to be degenerate, i.e. mi = m2. This choice is made for 
computational simplicity. 

III. OVERVIEW OF METHODOLOGY AND RESULTS 

In this section we discuss technical details relevant to the implementation of all the Ward 
identities, and give a summary of our results. 

The Ward identities can be implemented on states having any spatial momentum, and 
we collected data for (0, 1, a/2, \/3, 2) units of lattice momenta. In the extraction of ca 
we find that the results from all momenta are consistent, but only after errors proportional 
to {pay are taken into account. Because of these additional discretization errors, and the 
larger statistical errors in correlators with non-zero spatial momentum, we did not find that 
the results at non-zero momenta added useful information. Thus, for the calculation of all 
other renormalization and improvement coefficients we present results only from correlators 
with zero spatial momentum. 

We were unable to determine the covariance matrix to sufficient accuracy to do fully 
correlated fits. Thus, when fitting the time dependence of correlators, or ratios of correlators, 
we use only the diagonal part of the covariance matrix. Similarly, fits to the quark mass 
dependence (which are done within the jackknife procedure) ignore correlations between 
the results at different masses. Also, in the analysis of the three-point axial WI identities 
we do not propagate the errors associated with estimates of ca and cy as we do not have 
a corresponding error estimate on each jackknife sample. The fully self-consistent method 
would be to do a simultaneous fit to all the unknown parameters, but we do not have enough 
statistical power to do this. Because of these shortcomings, we can make no quantitative 
statement about goodness of fit. Nevertheless, assuming that the fits are good, the errors 
in the fit parameters, which are obtained using the jackknife procedure, should be reliable. 

In Table IH] we give our results for the critical hopping parameter k^, which is needed 
to define the vector Ward identity (VWI) quark mass m. These are obtained from the fits 
shown in Fig. ^ (for the 64NP dataset). We fit three quantities to quadratic functions of 
1/2k (which, up to an additive shift, is the tree-level quark mass). The first quantity fit is 
the quark mass m extracted from the axial Ward identity, using a mass dependent CA(m) 
and two-point discretization (see Section HVl below for definitions of these quantities). This 
is the middle curve in the plot. The second fit quantity is also fh, but now obtained using 
the chirally extrapolated value of ca- This is the lower curve in the plot. Finally, the third 
quantity fit is M^, and gives the upper curve in the plot. Only the last fit includes results for 
both degenerate and non-degenerate quarks, using the average value of I/k for the latter. 
As the curve shows, we find no noticeable dependence on the mass difference. The respective 
fit parameters are 



am{cA{rn)) 



am{cA{0)) 





(2/€)2 ■ 



(2k) 



1 



1 



2 ' 



(5) 



6 




1/2/c 

FIG. 1: Fits used to determine Kc by extrapolating 64NP results for m and M^. We show 
quadratic fits to the two-point version of fh for the two cases discussed in text (octagons label 
points with Cyi(m) and pluses label points with chirally extrapolated ca), and a quadratic fit to 
(diamonds). 

From these we get three estimates of Kc which we find to be consistent; this was not the case 
for 60NP and 62NP data. The first estimate, k''c\ is the most direct as fh and c^(m) are 
extracted together from the same two-point Ward identity (also see below), so we use it in 
subsequent analyses and, henceforth, drop the superscript. 

If both and fh are extracted from fits that include large times, where only the ground 
state survives, then it follows from eq. [7| below that 2m = M^{2/Bt, + qca) + O(a^), with 

oc (0|P|7r)//7r a quantity which is non-zero in the chiral limit (and which we will use in 
several places below). Thus M^- and fh should vanish at the same point. We use this fact 
to test the adequacy of our quadratic fits of versus fh{cA{fh)) or fh{cA{0)). The 64NP 
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60NP 






62NP 






64NP 




Label 




am 




K 


am 


aMj^ 


K 


am 


aMj^ 




0.1300 


0.1442(10) 


0.711(2) 


0.1310 


0.1345(6) 


0.609(1) 


0.1280 


0.2106(6) 


0.766(1) 


1^2 


0.1310 


0.1182(08) 


0.631(2) 


0.1321 


0.1053(4) 


0.522(1) 


0.1294 


0.1754(5) 


0.672(1) 


1^3 


0.1320 


0.0913(06) 


0.544(2) 


0.1333 


0.0728(3) 


0.418(1) 


0.1308 


0.1391(4) 


0.573(1) 


K4 


0.1326 


0.0752(05) 


0.488(2) 


0.1339 


0.0562(2) 


0.360(2) 


0.1324 


0.0960(2) 


0.449(1) 


K5 


0.1333 


0.0561(04) 


0.416(2) 


0.1344 


0.0419(2) 


0.307(2) 


0.1334 


0.0682(2) 


0.364(1) 


Kg 


0.1342 


0.0308(04) 


0.308(3) 


0.1348 


0.0306(2) 


0.261(2) 


0.1343 


0.0429(1) 


0.278(2) 




0.1345 


0.0236(35) 


0.265(12) 


0.1350 


0.0248(1) 


0.235(2) 


0.1348 


0.0285(1) 


0.223(2) 




0.13528(2) 






0.135854(5) 






0.135786(3) 






(2) 


0.13530(1) 






0.135875(4) 






0.135784(3) 






(3) 


0.13539(3) 






0.13594(2) 






0.13578(2) 







TABLE II: Values of the hopping parameter used in the various simulations, and the corresponding 
pseudoscalar mass aM^- and quark mass am defined using the Cyi(m) and two-point discretization 
(see Sec. II V|) . The three estimates of Kc, obtained using quadratic fits, correspond to (1) the zero 
of m with mass dependent ca_, (2) the zero of m with chirally extrapolated ca, and (3) the zero of 
Ml 



data, illustrated for two-point discretization of derivatives, give significant intercepts: 

a^M^ = 0.0090(23) + 2.61(5)am + 5.84(27)(am)2 (c^(m), /3 = 6.0) , 

a^M^ = 0.0098(27) + 2.54(7)am + 7.48(40)(am)2 (ca(0), /5 = 6.0) , 

a^M^ = 0.0049(13) + 1.87(3)am + 6.32(15)(am)2 (cA(m), /3 = 6.2) , 

a^M^ = 0.0043(13) + 1.86(3)am + 7.24(17)(am)2 (ca(0), /5 = 6.2) , (6) 

a^M^ = 0.0020(12) + L49(2)am + 6. ll(09)(am)2 (cA(m), /5 = 6.4) , 

a'^M^ = 0.0059(12) + 1.34(2)am + 7. 73(12)(am)2 (ca(0), /3 = 6.4) . 

Using CA{m) leads to smaller intercepts, and because of this we use am{cA{m)) (rather 
than am(cA(0))) when making chiral extrapolations in the subsequent analyses. We note, 
however, that the intercept is not small when converted into physical units (~ (160 MeV)^), 
and does not show any significant decrease with a. In this context it is important to note 
that the range of the fits in physical units is different in the three cases and the lightest 
"pions" are heavy. The range of pion masses in the three cases are 550 — 1500, 680 — 1770 
and 850 — 2900 MeV respectively. Thus neglected contributions from chiral logarithms, 
which become significant only at lower quark masses, and higher order terms in the chiral 
expansion, could account for the intercept. Since the present data are well fit by a quadratic, 
we cannot empirically resolve the issue of what additional terms need to be included in the 
fits. 

An important point to keep in mind is that the extrapolations to extract renormalization 
and improvement constants are in am, and are different from the usual chiral extrapolations 
where the control parameter is M^/A^ with A^ ~ 1 GeV. We do not need to be in the chiral 
regime for our method to work. In the ratios of correlators that appear in the Ward identities 
we use, the same intermediate states contribute to both numerator and denominator, and 
possible non-analytic behavior in the quark mass (including that from enhanced quenched 
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chiral logarithms) cancels. What matters is that am <^ 1, which, as can be seen from 
Table |nl is reasonably well satisfied for all of our masses. Indeed, a striking feature of our 
results is that the quadratic fits we use work very well over our entire range of quark masses. 





60NPf 


60NPb 


62NP 


64NP 




-0.039(08) 


-0.037(09) 


-0.034(03) 


-0.032(03) 


Zy 

bv 
bv 


+0.7689(08) 
+1.448(20) 
+0.7689(08) 
+1.530(12) 


+0.7703(09) 
+1.413(23) 
+0.7697(09) 
+1.519(13) 


+0.7880(04) 
+1.273(10) 
+0.7876(03) 
+1.402(08) 


+0.8033(05) 
+1.212(11) 
+0.8016(05) 
+1.370(09) 


_ '^V_ 

bA - bv 

Z'v 

bA - bv 

bA - bv 

PI Af'S 
bp - bs 

CT 


+0.773(05)* 
-0.309(76)* 
+0.774(05)* 
-0.288(69)* 
+1.203(11)* 
-0.094(89)* 
+0.799(04)* 
+1.052(10)* 
-0.058(65)* 
+0.083(12)* 


+0.761(06)* 
-0.469(77)* 
+0.762(06)* 
-0.433(71)* 
+1.209(13)* 
+0.016(96)* 
+0.794(05)* 
+1.062(12)* 
-0.178(52)* 
+0.088(12)* 


+0.790(02)* 
-0.096(31)* 
+0.791(02)* 
-0.095(29)* 
+1.191(06)* 
-0.075(75)* 
+0.818(04) 
+1.084(05)* 
-0.096(25)* 
+0.063(10)* 


+0.801(02) 
-0.123(54) 
+0.800(02) 
-0.130(49) 
+1.173(05) 
-0.131(65) 
+0.825(02) 
+1.089(04) 
-0.104(56) 
+0.054(05)* 


J%/Z\Zl[cA{m)] 
bA-bp- bm [cAim)] 

-2bm [cAim]] 
J^p/Z^ZllcAm 
bA-bp - bm [ca(0)] 
-2bm [ca(0)] 


+1.051(10)* 
+0.598(43)* 
+1.052(72)* 
+1.055(09)* 
+0.943(114)* 
+1.406(107)* 


+1.057(12)* 
+0.629(49)* 
+1.251(81)* 
+1.060(10)* 
+0.932(137)* 
+1.472(124)* 


+1.084(05)* 
+0.674(21)* 
+1.313(27)* 
+1.090(05)* 
+0.957(27)* 
+1.428(28)* 


+1.077(02) 
+0.511(09) 
+1.193(15) 
+1.077(03) 
+0.646(28) 
+1.346(13) 



TABLE III: Summary of results for the different combinations of renormalization and improvement 
constants extracted using the two-point derivative. The horizontal lines separate the extraction of 
quantities using the divergence of axial current, the conservation of charge, three-point axial chiral 
Ward identities, and the relation between quark masses given in Eq. 1201 All unmarked estimates 
are based on quadratic fits in both rhi = m2 and fh^. Asterisks mark values extracted using linear 
extrapolations in both mi = m2 and ms. All seven masses are used at /3 = 6.2 and 6.4 while at 
/? = 6.0 the lightest quark (^7) is dropped. Labels CA{rn) and ca(0) refer, respectively, to whether 
the mass dependent or chirally extrapolated value of ca is used in the analysis. 



With am and Kc in hand we carry out the analysis for two-point and three-point Ward 
identities discussed in Ref. [HI Each identity allows us to extract one or more combinations 
of on-shell improvement and normalization constants. Since many of the results for 60NP 
and 62NP data sets given in have changed as a result of our reanalysis, estimates from 
all three lattice spacings are given in Tables UTTl and HVl Similarly, a detailed comparison of 
the results for cy obtained using the methods discussed in Ref. |5| is given in Table IVl for all 
three values of 13. 
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60NPf 


60NPb 


62NP 


64NP 




— n (T^fih f\\ 

U.UOUl ±U ) 


u.uooi ±0 ) 






Zy 
bv 

Z'v 
bv 


+0.7677(25) 
+1.492(68) 
+0.7672(21) 
+1.533(13) 


+0.7679(31) 
+1.499(83) 
+0.7678(23) 
+1.522(14) 


+0.7877(04) 
+1.296(10) 
+0.7875(03) 
+1.402(08) 


+0.8027(06) 
+1.233(11) 
+0.8016(05) 
+1.370(09) 


bA - bv 

Z'v 

bA - bv 

bA — bv 

^A 
7O /7O 7O 

PI Af'S 

bp - bs 

CT 


+0.771(16)® 
-0.234(61)@ 
+0.773(16)® 
-0.210(56)® 
+1.198(08)® 
— u.u iy\t 'ij^ 
+0.800(04)® 
+1.051(08)® 
+0.003(42)® 
+0.087(11)® 


+0.760(17)® 
-0.359(57)® 
+0.762(17)® 
-0.324(54)® 
+1.196(08)® 

— \).yj\jz,\oo )^ 

+0.797(04)® 
+1.062(10)® 
-0.142(36)® 
+0.085(12)® 


+0.785(03)* 
-0.053(30)* 
+0.785(03)* 
-0.053(29)* 
+1.189(05)* 

— u.UO ( J* 

+0.814(02)* 
+1.084(05)* 
-0.096(25)* 
+0.071(10)* 


+0.797(02) 
-0.174(56) 
+0.797(02) 
-0.168(51) 
+1.173(04) 
— u.uuy ^^00 j 
+0.824(02) 
+1.087(04) 
-0.033(55) 
+0.058(05)* 


J%/Z\Zl[cA{m)] 

OA Op — Dm [CAyiil'Jl 

-2bm [cAim)] 
J^p/Z^ZllcAm 
bA-bp - bm [ca(0)] 

-2bm [ca(0)] 


+1.011(27)* 

+ 1.120(104)* 
+ 1.014(23)* 
-0.208(253)* 
+0.917(125)* 


+1.009(26)* 

+ 1.168(92)* 
+ 1.010(26)* 
-0.209(291)* 
+0.953(135)* 


+1.071(05)* 

+ 1.196(34)* 
+1.068(05)* 
+0.187(27)* 
+1.150(26)* 


+1.074(03) 

+1.183(19) 
+1.072(03) 
+0.233(31) 
+0.918(09) 



TABLE IV: Summary of results for the different combinations of renormalization and improve- 
ment constants extracted using the three-point derivative. All unmarked estimates are based on 
quadratic fits in both mi = rh2 and rhj,. Asterisks mark values extracted using linear extrapo- 
lations in both mi = m2 and ms. At /3 = 6.2 and 6.4, all seven masses are used. At /3 = 6.0 
the lightest quark K7 is dropped, and, for estimates marked @, only masses ki — K5 and linear fits 
are used. Labels CA{m) and Cyi(O) refer, respectively, to whether the mass dependent or chirally 
extrapolated value of ca is used in the analysis. 



Our final results for the individual constants are collected in Table IVII We quote both 
a statistical error (given by the single elimination jackknife procedure, in which we repeat 
the entire analysis on each jackknife sample), and an estimate of the residual 0{a) uncer- 
tainty. The latter is taken to be the difference in results obtained using two- and three-point 
discretizations of the derivatives except for hy where there is no three-point estimate. A 
different estimate of the uncertainties can be obtained by comparing our results to previous 
estimates by the ALPHA collaboration ^ ^ summarized in Table IVII by the QCDSF 
collaboration given in Table HXI Q and by the SPQcdR collaboration [15]. 

We collect separately, in Table IVni our results for the improvement constants c^, the 
coefficients of the equation-of-motion operators. These are discussed in Sec. IXllI In sec- 
tions IXIIII and IXIVI we present an analysis of residual discretization errors by comparing 
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64NP 




2pt 


3pt 


CAim) 


ca(0) 


CA{m) 


ca{0) 


extrapolation 


-0.167(90) 


-0.136(90) 


-0.042(71) 


-0.092(69) 


1/m fit 


-0.079(07) 


-0.090(08) 


-0.099(08) 


-0.079(06) 


slope ratio 


-0.088(10) 


-0.085(10) 


-0.066(08) 


-0.074(08) 


62NP 




2pt 


3pt 


CAim) 


ca(0) 


CA{m) 


ca(0) 


extrapolation 


-0.143(166) 


-0.128(164) 


+0.008(160) 


-0.036(157) 


1/m fit 


-0.104(17) 


-0.124(20) 


-0.161(20) 


-0.120(19) 


slope ratio 


-0.116(23) 


-0.117(23) 


-0.103(21) 


-0.105(21) 


60NPf 




2pt 


3pt 


CA{m) 


ca(0) 


CA{m) 


ca(0) 


extrapolation 


-0.058(157) 


-0.075(170) 


+0.158(115) 


-0.015(188) 


1/m fit 


-0.138(25) 


-0.189(47) 


-0.230(43) 


-0.143(58) 


slope ratio 


-0.136(23) 


-0.143(28) 


-0.092(16) 


-0.102(16) 


60NPb 




2pt 


3pt 


CA{m) 


ca(0) 


CA{m) 


ca(0) 


extrapolation 


+0.068(162) 


+0.112(186) 


-0.087(156) 


-0.036(211) 


1/m fit 


-0.132(23) 


-0.176(44) 


-0.221(43) 


-0.139(60) 


slope ratio 


-0.123(27) 


-0.128(35) 


-0.089(25) 


-0.104(26) 



TABLE V: Results for cy using the two-point discretization data. See text (sec lVH]) and Ref. [sf for 
details. The labels CA{m) and ca{0) refer to whether the mass-dependent or chirally extrapolated 
value of CA was used in the analysis. 



our estimates with those by the ALPHA collaboration and with one-loop tadpole improved 
perturbation theory estimates summarized in Table IVll 



IV. CALCULATION OF ca 

The calculation of ca exploits the two-point axial Ward identity 

EAdM, + acAd,P]^^^Kx,t)j(^HO)) _ 

which also defines the quark mass fhij. Here the superscript {ij) refer to the mass (flavor) 
labels {KiKj) of the quark and the antiquark. Up to corrections of O(a^), this ratio of 
correlators should be independent of the source J and the time t provided csw (the coefficient 
of the Sheikholeslami-Wohlert term in the action) and ca are tuned to their non-perturbative 
values. Since this criterion is automatically satisfied when the correlators are saturated by 
a single state, the determination of ca relies on the contribution of excited states at small t. 
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(3 


= 6.0 




/3 = 6.2 


/3 = 6.4 




LANL 


ALPHA 


P. Th. 


LANL 


ALPHA 


P. Th. 


LANL 


ALPHA 


P. Th. 


csw 


1.769 


1.769 


1.521 


1.614 


1.614 


1.481 


1.526 


1.526 


1.449 




+U.7Dy5(Sj(i9) 


+U.7»U9(d) 


+0.810 


-f-U. 1 O 1 o\^)\y^ ) 




1 n S91 




_uo so'^9^'f^^^'1 9^ 


_i_n s'io 




_i_n 7Q7('4V2'l 


_i_n 7Qnfi('94") 


+0.829 


+0.815(2)(1) 


+0.807(8)(2) 


+0.839 


+0.825(2)(1) 


+0.827(8)(1) 


+0.847 




+0.840(7)(15) 


+0.840(8) 


+0.956 


+0.883(4)(6) 


+0.886(9) 


+0.959 


+0.894(3) (4) 


+0.908(9) 


+0.962 




— U.UoDl^O )\0) 


— U.UOol^O ) 


— U.UiO 


— U.Uo4(oj(Dj 


— O.Uo8(4j 


—0.012 




— U.02o(2 ) 


—0.011 




— \J.\-O\A)\0 ) 


— U.oZl^ 1 ) 


— U.Uzo 




—0.21(7) 


—0.026 


-0.09(1)(2) 


—0.13(5) 


—0.024 






NT A 


in 090 




IN .A. 


+0.019 


H-u.Uo4(oj(4j 


IN .A. 


+0.018 


bv 


+1.43(2)(7) 


N.A. 


+1.106 


4-1.27(l)(2) 


N.A. 


+1.099 


+1.21(1)(2) 


N.A. 


+1.093 


bv 


+1.52(1) 


+1.48(2) 


+1.274 


+1.40(1) 


+1.41(2) 


+1.255 


+1.37(1) 


+1.36(3) 


+1.239 


bA - by 


— 0.21(6){5) 


N.A. 


-0.002 


-0.09(4)(2) 

\ /V / 


N.A. 


-0.002 


-0.13(5)(6) 


N.A. 


-0.002 


bA — by 


-0.36(7){13) 


N.A. 


-0.002 


-0.10(3)(8) 


N.A. 


-0.002 


-0.13(5)(3) 


N.A. 


-0.002 


bp — bs 


-0.12(6){5) 


N.A. 


-0.066 


-0.10(3)(1) 


N.A. 


-0.062 


-0.10(6)(7) 


N.A. 


-0.059 


Ia - bp +bs/2 


+0.61(5)(12) 


N.A. 


+0.585 


+0.67(2)(6) 


N.A. 


+0.579 


+0.51(1)(1) 


N.A. 


+0.575 


bs 


+1.15(8)(1) 


N.A. 


+1.172 


+1.31(3)(12) 


N.A. 


+1.161 


+1.19(2)(1) 


N.A. 


+1.151 


bA 


+1.22(6){11) 


N.A. 


+1.104 


+1.19(4)(5) 


N.A. 


+1.097 


+1.09(5)(6) 


N.A. 


+1.092 


bA 


+1.16(7)(10) 


N.A. 


+1.271 


+1.31(3)(4) 


N.A. 


+1.252 


+1.24(5)(4) 


N.A. 


+1.237 


bp 


+1.02(9)(15) 


N.A. 


+1.105 


+1.19(4)(3) 


N.A. 


+1.099 


+1.13(6)(1) 


N.A. 


+1.093 



TABLE VI: Final results for improvement and renormalization constants. The first error in LANL 
estimates (this work) is statistical, and the second, where present, corresponds to the difference 
between using 2-point and 3-point discretization of the derivative. We quote both by, bA and by, bA 
to simplify comparison with previous results. Estimates for Zp/Zg by the ALPHA collaboration 
are taken from Ref. and the rest from Refs. [ill . H^ . [3- The final column gives the results 
from one- loop tadpole improved perturbation theory Ref. [5|. 



The sensitivity of the ratio ((Zj) to ca is illustrated for the 64NP data in Fig|2l We find 
that, for J = P, the contribution of higher excited states is significant only at time-slices 
t = 1 — 5 for all three values of /3 (see p| for data at P = 6.0 and 6.2). For two-point and 
four-point discretization, the data at t = 1 cannot be used to extract ca as the discretization 
of d\P{t = 1) in Eq. [7| overlaps with the source at time slice t = 0. Consequently, only the 
range 2 < t < 5 is sensitive to tuning ca and we choose ca to make fhij as flat as possible 
for timeslices t > 2. This is done by minimizing the for a fit to a constant, as illustrated 
in Fig. 121 For three-point discretization we can implement the same choice only for t > 3. 
The J = A4 data are not presented as they are dominated by the ground state already at 
t < 4 and are thus not sensitive to the choice of ca- Results for m with different choices of 
discretization are collected in Table IVlTTl 

As noted in Ref. a possible problem with our criterion for determining ca is that 
it does not involve the same physical distances at all couplings. The "physical" criterion 
we wish to implement is that the same value of fhij is obtained in the AWI from both the 
ground and the first excited states. This requires that we tune the source to produce the 
same mixture of ground and excited states at all values of a. While the lattice size and the 
radius of the smeared source in the generation of quark propagators were increased with /3, 
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60NPf 


60NPb 


62NP 


64NP 


dy + Cp 


2.97(26) 


2.66(36) 


+2.68(9) 


2.51(13) 


+ Cp 


2.56(20) 


2.54(24) 


+2.44(10) 


2.23(12) 


2c'p 


2.98(41) 


2.76(43) 


+2.99(13) 


2.12(22) 


c'^ + Cp 


2.54(14) 


2.54(15) 


+2.38(7) 


2.36(11) 




2.58(19) 


2.48(23) 


+2.45(11) 


2.40(13) 


c'v 


1.48(31) 


1.28(32) 


+1.19(10) 


1.45(13) 


c'a 


1.07(27) 


1.16(24) 


+0.94(11) 


1.17(11) 


C'p 


1.49(21) 


1.38(22) 


+1.50(7) 


1.06(11) 


c's 


1.05(22) 


1.16(19) 


+0.89(9) 


1.30(10) 




1.09(25) 


1.10(22) 


+0.95(12) 


1.34(12) 



TABLE VII: Results for off-shell mixing coefficients using the two-point derivative data and Cyi(m). 







60NP 






62NP 






64NP 




Label 


2-pt 


3-pt 


4-pt 


2-pt 


3-pt 


4-pt 


2-pt 


3-pt 


4-pt 


(k1, Kl) 

(k2, K2) 

(«;3, /^s) 
(k4,K4) 


0.1442(10) 
0.1182(8) 
0.0913(6) 
0.0752(5) 


0.1425(13) 
0.1167(11) 
0.0902(8) 
0.0744(7) 


0.1450(10) 
0.1188(8) 
0.0920(6) 
0.0757(5) 


0.1345(6) 
0.1053(4) 
0.0728(3) 
0.0562(2) 


0.1340(7) 
0.1048(5) 
0.0722(3) 
0.0558(3) 


0.1349(6) 
0.1058(4) 
0.0727(3) 
0.0561(2) 


0.2106(6) 
0.1754(5) 
0.1391(4) 
0.0960(2) 


0.2108(9) 
0.1752(7) 
0.1389(5) 
0.0958(3) 


0.2110(6) 
0.1756(5) 
0.1394(4) 
0.0963(2) 


(^5,^5) 
(^6,^6) 
(k7,K7) 


0.0561(4) 
0.0308(4) 
0.0236(35) 


0.0555(6) 
0.0312(15) 
0.0251(55) 


0.0563(4) 
0.0310(4) 
0.0246(30) 


0.0419(2) 
0.0306(2) 
0.0248(1) 


0.0417(2) 
0.0304(2) 
0.0247(2) 


0.0421(2) 
0.0307(2) 
0.0250(1) 


0.0682(2) 
0.0429(1) 
0.0285(1) 


0.0681(2) 
0.0428(1) 
0.0284(1) 


0.0686(2) 
0.0430(1) 
0.0286(1) 



TABLE VIII: Values of the quark mass am defined using the mass-dependent ca for two-point, 
three-point and four-point discretization schemes at the three couplings. 

they were not tuned. In fact, we find that our fit is sensitive to the same range, t = 2 — 5 
(3 — 5 for three-point discretization), for all three lattice spacings, and thus is sensitive to 
significantly shorter Euclidean times at /3 = 6.4 than at /? = 6. 

The part of our analysis which is, therefore, sensitive to the extent to which our criterion is 
physical is the manner in which the continuum limit is approached. With a physical criterion, 
the dominant correction to scaling will be proportional to 0?. If one changes the criterion 
as a is varied, the simple power dependence can be distorted. While the data suggests that 
the contribution from higher states is small, we cannot rule out some distortion and the 
scaling analysis has to be taken with caution. To study the question in detail, however, 
would require a more extensive data set than ours. 

In Figs 01 m and El we show quadratic fits to Cyi(m) versus rh for the two-point, three-point 
and four-point discretization of the derivative. The data are for zero momentum correlators 
at /9 = 6.0, 6.2 and 6.4, and include degenerate and non-degenerate mass combinations. 
The results for two-point and three-point discretizations are given in Tables 11111 and IIVI 
The four-point estimates are 0.034(8), 0.034(9), 0.030(4) and 0.030(3) for 60NPf, 60NPb, 
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62NP and 64NP data sets respectively. We find that the three estimates agree within 
errors in the chiral hmit. By contrast, the 0{am) contributions are significant, as shown by 
the large, roughly linear, dependence of ca on fh. However, as we now explain, the bulk of 
this linear dependence has a simple kinematic origin and can be understood analytically. 

In Ref. we showed that if, by tuning ca and m, eq. ((Tj) can be satisfied over a common 
range of time-slices where two- and three-point discretizations schemes are implemented 
then, to O(a^), in is the same in both schemes and the ca at any am are related as c^~^* = 
c^~^* — ma/2 + (a^). It is useful to generalize this argument to provide the relation between 
Ca determined in any two discretization schemes. The equation we wish to satisfy is (94^4 + 
acAd\P—2mP = 0. Taylor expansion of any lattice version of this relation using a symmetric 
discretization scheme for the derivatives gives 

[d4A^ + a^a^^A^ + 0{a'^)] + acA[dlP + a^(3dlP + 0{a'^)]-2m[P + a^jdlP + 0{a'^)] = 0. (8) 

The key assumption is that this relation can be satisfied by two discretization schemes over 
a common range of timeslices. If so, then it follows, first, that the two schemes will give the 
same m to 0{a^) and, second, that 

cj) - c^i) = 2am(a(i) - a^^^ - 7^^) + 7(2)) + 0{a^) . (9) 

For the two-point and three-point derivatives we have used, this condition reduces to c^~^* = 
4"^* - am/2 + 0(a2) because a^-p* = 1/24, 72-P* = 1/8, a^-P* = i/g, and 7^^^* = 0. Our 
data confirm these two predictions to good accuracy: the ratio of correlators (2mjj) is the 
same within errors for the three discretization schemes, as shown in Table IVIIIl and Eq. |H1 
holds, as illustrated in Figs. El IH and El where we also plot the quantity c^~^* — c^~^* + am/2 
and show it is consistent with zero at the 1-a level. 

The UKQCD collaboration has pointed out that the mass dependence of ca can be 
reduced using higher order discretization schemes. This also follows from Eq. El For any 
O(a^) improved scheme, = 7 = 0, and consequently = c'^^^ — am/Q + O(a^). In fact 
we find that the slope of c^~^* versus arhij is A ^ 0.18, 0.18 and 0.19 respectively for /3 = 6, 
6.2 and 6.4, so that the slope obtained using any O(a^) improved scheme should indeed 
be very small, A — 1/6 ~ 0.02. For our four-point discretization scheme (which is O(a^) 
improved, but differs from the five-point scheme used in Ref. 0) we find that the slope is 
< 0.03 for all four data sets as shown in Figs El IH and El We also find that the contribution 
of the {amf' term in the four-point scheme is comparable to the linear term over the range 
of quark masses simulated, and to the errors. Because of the size of these higher order terms, 
further improvements in the discretization of the derivative are not expected to reduce the 
undertainty. 

We stress that we do not expect a higher order scheme to completely remove am contri- 
butions in cai for to do so would require complete implementation of an O(a^) improvement 
program.^ Nevertheless, our results indicate that the bulk of the slope for two- and three- 
point discretization is due to errors associated with discretization of the derivative. 

The upshot of this discussion is as follows. On the one hand, it would have been ad- 
vantageous to use a higher-order discretization scheme with a smaller slope A. This would 



^ This means that the assumption leading to eq. ©, namely that the relation JHJ can be satisfied by two 
schemes over a range of timeslices, cannot hold precisely. 
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have reduced the uncertainty in the results for some of the bo that are proportional to A, 
as discussed in later sections. On the other hand the 0(am)^ and 0{am) terms become 
comparable in our four-point data, and the error in the extrapolated value does not decrease 
compared to the lower-order schemes.^ Having demonstrated that the dominant effect is 
kinematic, we can remove it in our two- and three-point schemes simply by using the mass- 
dependent c^(m) in the improved axial current, rather than Ca(0). This is indeed what we 
do in the axial rotation 5S (which, we recall, is always defined with two-point derivatives 
and in which we always use the two-point ca(w^)). 

We now return to the numerical results for ca(0), which are given in Tables IIIII and 
IIVI Results from two-, three- and four-point discretizations should differ only by terms of 
relative size aAqcn- In fact, as already noted. Figs 01 ID and El show that the quadratically 
extrapolated values for ca at each of the three jS from all three discretization schemes agree. 
Estimates using linear fits, however, differ by combined la errors, due to the curvature. 

In Ref . i we chose, for our central value, ca from the two-point discretization of the 
derivative over that from the three-point derivative for the following two reasons. First, the 
O(a^) discretization errors in the derivatives are smaller, which leads to a smaller slope of 
Ca versus m; and second, because the statistical errors are smaller. Now that we understand 
the relative size of the slope to be largely a kinematical effect, and the extrapolated values 
overlap, and the uncertainty in the estimates are comparable, we take the weighted mean 
of the two-, three- and four-point results for our central value at all three lattice spacings. 
In addition to statistical errors we quote the spread of the results to estimate the residual 
0(a) errors. Note that, as pointed out in the Introduction, the choice of the discretization 
scheme used for the derivative does not affect results for c^(0) at the leading order of 
overall improvement, and the estimates from any scheme can be used to define the improved 
theory. However, if a calculation requires the axial vector Ward identity be respected, then 
the appropriate discretization scheme and the corresponding CA{m) should be used. 

V. EXTRACTING ca USING STATES AT FINITE MOMENTUM 

The data at f3 = 6.2 and 6.4 are precise enough to extract ca using states having non-zero 
momenta. In figures (HI and |7| we show the results of linear fits for the chirally extrapolated 
value of Ca for two-, three-, and four-point discretization of the derivative as a function 
of (pa)^. For the chiral extrapolation of ca(^), quadratic fits to all mass combinations of 
Ki — K,7 propagators work very well at all five values of pa. The signal in three-point data 
at pa = 2 is noisy and this is reflected in the errors. The data exhibit the following two 
features: 

• Additional discretization errors of O(p^a^) are generated when using states with non- 
zero momenta. The coefficients of these corrections are significant, lying in the range 
0.12 - 0.22. 

• The difference in results between the two-, three- and four-point discretization of the 
derivative decreases significantly between j3 = 6.2 and 6.4. 



^ Unfortunately, we are not able to determine the efficacy of the four-point discretization scheme for analyz- 
ing the three-point Ward identities as some of the required raw data has been lost due to disk corruption. 
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FIG. 2: Estimates of 2afhij for different values of ca illustrated using i = j = k^, J = P, and 
two-point discretization in the 64NP data set. For this quark mass, ca = —0.0089 extends the 
plateau to the earliest allowed time slice t = 2. To show sensitivity to the tuning we contrast this 
best fit with those using ca = and ca = —0.03, the latter being close to the chirally extrapolated 
value. 



Overall, the consistency of the estimates between the two-, three- and four-point discretiza- 
tion schemes, the added information from fits versus (pa)^, and the expected improvement 
with P enhance our confidence in our quoted estimate of ca- 



VI. Z'^ AND bv 

Our best estimate of Zy comes from the matrix elements of the vector charge 
J d^xV}'^^\x) between pseudoscalar mesons 

1 _ E,,y{P^''K^,r){Vj)f\y,t)j(^^\0))_ ^^^^ 



with r > t > and the superscript (23) denoting the flavor of the two fermions in the bilinear 
(which are taken to be degenerate). Our results for this ratio, illustrated by those in Fig. |H1 
show two features of particular interest: first, there is a significant dependence on t for time 
slices close to the source or sink; and, second, there is a clear difference between the results 
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FIG. 3: Results of chiral extrapolation of ca data at /3 = 6.0 (60NPf) for the two-point, three- 
point and four-point discretization of the derivative. Quadratic fits are made to all degenerate and 
non-degenerate mass combinations using ki — kj. We also show the quantity c^~^* — c^~^* + afh/2 
discussed in the text using the symbol plus. 

using J = P and J = A4. Both features are indicative of O(a^) corrections, since, aside from 
such corrections the ratio should be independent both of t and the choice of source. The 
observed effect is (1 — 5) o^Aqqjj using Aqcd = 300 MeV. Since neither feature would be 
present if the source and the sink coupled to a single state (irrespective of improvement), our 
results show that the separation between source and sink in our calculation is insufficient to 
isolate the lowest state for any value of t. We stress, however, that this is not a problem for 
implementing the improvement program (since, after all, we expect ambiguities of 0(a^)). 

To obtain our central values we average the J = P and J = A4 data within the jackknife 
procedure as they are of similar quality. There is a slight difference in results for Zy using 
two- and three-point derivatives, as shown in Tables UTTl and HVl Also, ai (3 = 6.0, the errors 
in the three-point estimates are almost three times as large. These differences arise during 
the chiral extrapolation because the rh and Kc are slightly different for the two cases and 
the errors in c^(m) are 2 — 3 times larger for the three-point data. The difference between 
linear and quadratic chiral extrapolation, as shown in Fig. |H1 is significant. For our central 
values we use quadratic extrapolations of the two-point data. 

To extract Zy, by and by we fit the ratio in eq. (|10|) to a quadratic function of both fh 
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FIG. 4: Same as Fig. El but at /3 = 6.2. 

and the VWI mass. At (3 = 6A the fits yield 

Zv = 0.8033(5) [1 + 1.212(ll)ma + 1.134(39)(ma)2] , (11) 
Zv = 0.8016(5) [l + 1.370(9)ma + 0.033(24)(ma)2] . (12) 

The two intercepts, which give Zy, differ by 3-a at /? = 6.4 and by < 1-cr at /? = 6.0 and 6.2. 
For the final estimate of Zy we choose the weighted average as they are of similar quality. 
The coefficient of the linear term in the two fits gives by and by respectively. 

In Table IIXI we compare results with those from other non-perturbative calculations 
that have been done with the same 0{a) improved fermion action but utilizing different 
initial and final states to measure the charge. We find that the results for by agree within 
the combined la uncertainties and the expected differences of 0(a). For Zy, there is a 
significant difference between the LANL and the ALPHA 0, 0, ^| collaboration values, 
which we show, in section IXIIH can be explained as residual O(a^) effects. Estimates by the 
QCDSF [13] and the SPQcdR \xE] collaboration lie in the range defined by the LANL and 
ALPHA data. 

In Ref. jif it was observed that extrapolations using a quadratic fit in m give estimates 
closer to measured values of Zy near the charm quark mass than did fits versus m. At 
f3 = 6.4 the two fits agree within 1% up to ma ^ 0.3, whereas the charm quark mass is 
smaller in lattice units, i.e. arric ~ 0.22. Since, as noted above, O(a^) errors are ~ 1%, we 
conclude that either fit can be used for quark masses in the range — rric. 
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FIG. 5: Same as Fig. El but at /3 = 6.4. 





(3 = 6.0 


(3 = 6.2 


/? = 6.4 




LANL 


ALPHA 


QCDSF 


LANL 


ALPHA 


QCDSF 


LANL 


ALPHA 


QCDSF 




0.7695(8) 


0.7809(6) 


0.7799(7) 


0.7878(4) 


0.7922(4)(9) 


0.7907(3) 


0.8024(5) 


0.8032(6)(12) 


0.8027(2) 


bv 


1.52(1) 


1.48(2) 


1.497(13) 


1.40(1) 


1.41(2) 


1.436(8) 


1.37(1) 


1.36(3) 


1.391(5) 



TABLE IX: Non-perturbative estimates of and by from the LANL, ALPHA, and QCDSF Q 
collaborations. For consistency LANL estimates at all three (3 are taken from fits versus the VWI 
mass m. 

VII. cy AND 6a -V 

Up to this stage, we have used only two-point correlation functions or three-point correla- 
tors involving the vector charge. We now turn to axial Ward identities involving three-point 
correlators. These allow us to determine cy, Ba — by, Z^, Zp/Z^, bp — 65, and ct, as well as 
giving an alternate determination of Zy. We first consider the improvement coefficient cy 
whose precise determination feeds into the calculation of Z^, Zp/Zg, ct, and c^- The best 
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FIG. 6: CA using states of non-zero momentum plotted against (pa)^, along with a linear fit. We 
show results for the two-point, three-point and four-point discretization of the derivative. The data 
are for P = 6.2 and pa is in units of 27r/24. 



signal for cy is obtained by enforcing Ni = N2 + cyD, with 



so that 



D 



(23), 



l(31), 



CV 



N 
D 



AWy2 

D 



(13) 
(14) 
(15) 

(16) 



We recall that 6S uses two-point discretization (and the corresponding value of Cyi(^i)), but 
that the other improved currents in these expressions are discretized both with the two- and 
three-point forms giving two sets of estimates. Within each set we provide two estimates 
using ca{0) and CAifn) in the expression for Aj. We also recall that we always use rhi = m2. 
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FIG. 7: CA using states of non-zero momentum at /? = 6.4. Notation as in Fig. IHl 

Figures ITIH ITD and IT2l illustrate the quality of our data for A^^i, N2 and D. The improve- 
ment in errors and overall quality as jS increases is evident. Note that N2 and D are expected 
to have larger errors than Ni since the lightest state which contributes is the axial-vector 
rather than the pion. 

Our procedure is to determine Ni, N2 and D from fits to the plateaus and then combine 
the first two to form N = Ni — N2. The results for N and D at (3 = 6.4 are shown in Fig. IT^ 
where it is apparent that the errors in N determine the quality of the result for cy = N/D. 
As noted in Ref. @ for the data at /3 = 6 and 6.2, both and D are to good approximation 
functions of rhi — that vanish when mi ^ m^. Since they do not, however, vanish at 
exactly the same point (presumably due to statistical and residual discretization errors), 
their ratio diverg shown in Fig. El 

In Ref. i, we used three methods to extract cy that try to minimize the effect of this 
spurious singularity, and we follow the same strategy here. Details of the methods will not, 
however, be repeated. Our estimates at /3 = 6.0 and 6.2 have changed after redoing the 
chiral fits to A^i, A^2 and D, and so we quote, in Table IVi results for all (3. For each method 
we have an additional four choices: we can use two-point or three-point discretization of the 
currents, and for each of these we can use either mass- dependent or chirally extrapolated 
values of ca in the operator {Ai)'^^^'^ appearing in the denominator of Ni. 

The extrapolation method has the largest uncertainty so we discard it. The consistency 
between the result using the "1/m fit", shown in Fig. 1141 and the "slope-ratio" method. 
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FIG. 8: The ratio defined in eq. Ellfor sources J = P and J = A4. Data from the 64NP set with 
ah propagators having mass k = K5. 

improves with /3, but the "slope-ratio" method is more stable with respect to the range of 
quark masses used in the fits at all three (3 values, and has the smallest dependence on 
the choice of ca- We therefore take our final estimates from the "slope-ratio" method and 
average the Cyi(m) and CAi^h = 0) values to get our final estimates. As usual, we take the 
central value from the two-point scheme and use the three-point scheme to estimate the 
discretization error. 

We can also use the quantity A^i, defined in Eq. UHl to determine bA — by, by — bA and to 
give an alternate determination of Zy. We must first extrapolate to rhi = m2 = to remove 
the contribution of equations-of- motion operators. In Fig. [THlwe illustrate the quadratic fits 
used to do this for the 64NP data. We then fit to a quadratic function of fh^ or m^. These 
fits, shown in Fig. for two-point discretization and ca(0), have parameters 
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^(l + (6^-M^ + 0(a^)) = 1.249(3)(1 - 0.123(54)^ + 0.06(38)(^)^) 
^(l + (6^-M^ + 0(a^)) = 1.250(3)(1 - 0.130(49)^ - 0.25(33)(^)^) .(17) 

The estimates for Zy are consistent with those obtained using the conserved vector charge, 
Eq. ^2 but have larger errors, so our preferred value is from the analysis presented in 
Section ED 
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FIG. 9: Linear and quadratic fit to Zy versus am for the 64NP data set. 

The coefficient of the term hnear in m (m) gives 6a — 6y (6^ — fey). We find that the 
errors in both 6^ — &y and 6^ — &y are large and comparable. In addition, there can be large 
0(a) errors feeding in from the dependence of ca on m as discussed below. 

It is easy to see that when using Eq. (fT^ to extract hA — by the result will depend on the 
choice whether Cyi(m) or the chirally extrapolated ca{0) is used. As explained in a shift 
Ca ^ Ca + ^ in the definition of (^4^)4 in the denominator produces a change in Zy of the 
form ^aBj,. If, instead, we use CA{fh) = ca + Ama in the calculation then the slope, not the 
intercept, changes, i.e., one gets feyi — 6y — Aai?^/2 instead of 6^ — fey. For the two-point data 
at /? = 6.4 the two estimates are —0.32(5) and —0.12(5) for CA{m) and ca(0) respectively. 
The difference, ~ 0.20, even though formally of higher order in a, is large because aB-j^ ~ 1.5 
and A = 0.19 as discussed in the extraction of ca- We do not have an a priori argument 
that, to this order, favors one choice over another. Anticipating that calculations of physical 
quantities will use improvement constants defined in the chiral limit and understanding that 
the slope A is almost entirely an artifact of the discretization scheme used to calculate ca, 
we take results obtained using CA{rn = 0) for the two-point discretization as our estimates. 
We stress that we do not include the difference between the results using CA{rn) and ca{0) 
as part of the error. These new results supercede those given in j^. 

Overall, Ba — by is small and the uncertainty is comparable to the signal. The expected 
relation (pA -by) = {Z\Zl/Z%){bA - by) + 0(a) holds at the la level. 
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FIG. 10: Illustration of the quality of the signal for the quantity A'^i of Eq. (|13jl for all four data 
sets with two-point discretization and CA{fn). In all cases the data have to be multiplied by the 
respective values of 2^13, the lattice normalization of the additional propagator in the numerator. 



VIII. 

The Ward identity 

EyASSj'"^ {Ai)f'\y,y,)v!^'\0)) ^ Z'y{l + hyam,/2) 
Y.ymT\y,y,) v/^^)(0)) Zl . Z\{1 + hj,am,/2) ' 

gives Zy/{Z^'^ and a second estimate of 6a — by- The quality of the signal for the ratio of 
correlation functions, as illustrated in Fig. [TTl is good as the intermediate state is the vector 
meson. Data in Fig. show that quadratic fits in both rhi = m2 and are preferred 
aX (3 = 6.4. Linear fits are sufficient at /5 = 6.0 and 6.2. The resulting values are given in 
Tables Em and CVl 

Including the results in Section IVTl we have two estimates for — by with similar errors. 
These estimates come from Ward identities that involve different, pseudoscalar versus vector, 
intermediate states. Also, in Eq. ^Jthe term proportional to ca in Aj does not contribute at 
zero momentum so there is no associated uncertainty. Thus, the 0{a) errors can be different 
in the two cases. As shown in Tables IIIII and IIV[ we find that the two estimates show 
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FIG. 11: Illustration of the quality of the signal for of Eq. (|14j) for all four data sets with two- 
point discretization and Cyi(m). In all cases the data are for K3 and have to be multiplied by the 
respective values of 2^3, the lattice normalization of the additional propagator in the numerator. 



considerable 0{a) variation, but this is not unexpected given the size of the errors and the 
possibiUty of additional 0{aAQCD ~ 0.2 — 0.1) uncertainty in previous estimates as discussed 
in Section IVTl Had we chosen to use Cyi(m) to extract Ba — hy in Section IVll the variation 
would have been larger by a factor of two or more. Thus, for our final estimate we average 
the two two-point estimates and quote the difference between two-point and three-point 
discretization schemes as an estimate of residual 0{a) errors. The upshot of the analysis is 
that hA — by is small and the systematic errors are of the same size as the signal. 

To estimate we use the product of Eqs. (fTHj) and (fT^ as it yields 1/(Z°)^ directly. 
The final chiral extrapolation in arhs for the product is shown in Fig. ^| In this product the 
terms proportional to am^ cancel, but nevertheless the data show a clear ams dependence. 



This we interpret as due to O(a^) terms of the generic form Z^(m3) 



Z0(0)(l + a2m3A). 



The slopes at /3 = 6.0, 6.2 and 6.4 are 0.22(6), 0.12(11) and 0.11(4) respectively. To match 
the observed slope aA ^ 0.11 at /3 = 6.4 requires A ~ Aqcd ~ 0.4 GeV, which is a reasonable 
value. Also, the change between f3 = 6.0 and 6.4 is consistent with the expected scaling in 
a. In Ref. [5| we had ignored this dependence and fit the data to a constant to extract Z^. 
In light of our results at /3 = 6.4 and a better understanding of possible rhs dependence, we 
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FIG. 12: Illustration of the quality of the signal for D of Eq. (jlSf) for all four data sets with 
two-point discretization and c^ifn), using ^3 propagators in all cases. 



have refit the data at f3 = 6.0 and 6.2 also. We now use quadratic extrapolation in rhi = m2 
at /3 = 6.4 and linear dX (3 = 6.0 and 6.2. Linear extrapolation in works well at all three 
couplings, however at /3 = 6.2 and 6.4 we use quadratic fits to maintain consistency with the 
rest of the analysis. At these weaker couplings linear and quadratic estimates are consistent. 

A comparison of Figs. ^1 and ^1 raises the following concern. The slope in Fig. with 
respect to arh^ relative to the intercept is an 0(a) effect, proportional to hy — b^, while that 
in Fig.^Jis, as just discussed, of one higher order.'' The two slopes are, however, numerically 
very similar. This once again suggests that there can be substantial uncertainty of 0{a), 
comparable to the value itself, in any result for by — bA- In fact, our analysis illustrates 
a problem common to the extraction of all measurements of the differences ba — bso- The 
signal, the errors, and the O(a^mA) uncertainties are all comparable. 



^ Even though the data in Fig. Ells consistent with no ma dependence, this is not true at (3 = 6.0 and 6.2. 
We make a quadratic fit as indicated by all other data at /3 = 6.4. 
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FIG. 13: Data for N and D, defined in the text and used to extract cy, plotted as a function of 
nil — fnz for the 64NP dataset with two-point discretization and CAirh). 

IX. Z%/Z%bp-hs 

To obtain Zp/(Z^Z°) and hp — hs we use the identity 

evaluated in the limit m\ = m2 with J = P or ^4. The intermediate state in both the 
numerator and the denominator has the quantum numbers of a pion, and the ratio has a 
very good signal, whose quality, as a function of /3, is shown in Fig. 1201 Data at /? = 6.4 for 
the ratio on the left hand side of Eq. (fT^ favor quadratic fits for rhi = 1112 and — >■ 
extrapolations as illustrated in Figure EU The intercept and the slope give Z%/{Z^Z^j^ and 
hp — hs respectively, and these estimates are quoted in Tables ITTTl and HVl To get Z^/Z^ 
we eliminate Z\ by combining the ratio in Eq. (fT^ with the product of Ward identities 
discussed in section I Villi 

The value of hp — hs is numerically small, comparable to the errors and of the same 
order as 0{a?m^K) effects discussed in Section rVIIII In this case we take the average of the 
two-point and three-point values as our best estimate. The reason is that the operators in 
Eq. (|19|) do not contain any derivatives (no improvement terms) so the difference between 
two-point and three-point estimates arises solely from the chiral extrapolations due to the 
tiny differences in fh for the two cases as shown in Table IVIIII 

Our results for Zp/Zg, obtained using Eq. (fT^ and Za from Section IVIIH are presented 
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in Table IVII These are consistent with the recent estimates by the ALPHA and SPQcdR 
collaborations In Section fXIVI we compare our results with predictions of one- loop 

perturbation theory and discuss the size of O(a^) and 0{a'^) corrections needed to explain 
the large difference. 

X. ZO bA-bp- AND bm 

One can derive a relation between the two definitions of quark mass |1 

- = - (bA -bp -K)am,. + O^amf] , (20) 

m Z\ mav 

where Xav = {Xi+X2)/2. This relation is useful because = 1/Z° and bs = -2bm BQ. 

In Fig. 1^ we illustrate fits to Eq. (j^UI) for the simpler case of degenerate quarks for m 
calculated using both CA{m) and ca{0). In this case the term proportional to does not 
contribute. The data show that for (3 = 6.4 including a term quadratic in rhav gives a much 
better fit whereas for 6.0 and 6.2 a linear fit suffices. The term proportional to b^ contributes 
only to non-degenerate combinations. Fits using all combinations of six {(3 = 6.0) and seven 
(/3 = 6.2 and 6.4) values of quark masses allow both bp — bA — b^ and bm to be extracted 
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FIG. 15: Quadratic extrapolation of the ratio in Eq. (jlcij) in fhi = m2 for fixed fh^ = for 64NP 
data set with two-point discretization and Cyi(m). 

reliably. These two sets of results of the fits, using ca(^) and ca(0), are quoted in Tables ITTTl 
andEVl 

The intercept, which gives Zp/Z'^Zg, should be same for ca(?ti) and c^(0), to the extent 
that the fits are good. Furthermore, the difference between two- and three-point results 
(which have different results for ca) should be small. These features are borne out by the 
results. The only notable difference is that the errors in the two-point data are smaller. The 
results are also consistent with those obtained using Eq. (fT^ . and have similar errors, as 
illustrated in Fig. |^ In Ref. we preferred the results from Eq. (fT^ since the method 
of this section has a greater sensitivity to uncertainties in (which are enhanced by the 
presence of the factor B^^ = M^/m ^ 4 GeV). With better understanding of the errors we 
now choose to take, for our final value of Zp/Z^Z'^, the weighted mean of the two-point 
results from Eq. (fT^ and those using m/m, with the latter determined using CA{rn). 

The extraction of bA — bp + bs/2 and 65 is effected by the choice of ca- Using Eqs.rZlandEUl 
one can show that, to leading order, (&yi — &p + &s'/2)|c^(m) = (6a — &p + &5/2)|c^(o) — 
and similarly for bm- Our data are roughly consistent with this relation for both the two-point 
and three-point discretization methods. For example, in case of the two-point discretization 
method, the values of the slope A, illustrated in Fig. El for P = 6.4 data, are approximately 
0.18, 0.18 and 0.19 and ai?,^ ~ 2.6, 1.9 and 1.5 at the three couplings respectively. This 0{a) 
effect, enhanced by the large value of B^,, gives rise to the difference in slopes as illustrated 
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FIG. 16: Quadratic fits to Eq. (|13() to extract bA — by and bA — by for 64NP data set with two-point 
discretization and c^(?n). The crosses show data and fits versus the AWI quark mass m, whereas 
the diamonds show results obtained using the VWI quark mass m. 

in Fig. ini 

For the two combination of 6's the two-point and three-points results are consistent for 
c^(m). This is expected because, as discussed in section HVl at each quark mass the in 
extracted from the two discretization schemes are, up to O(a^), the same, provided the mass 
dependent c^(m) are used. We also find that the fits to two-point data with CAim) are 
marginally better. So we use estimates obtained from the two-point data with CAiin) for 
the central values. 

Note that the considerations regarding choice of CA{rn) versus ca{0) here are different from 
those applied in section FV 111 when determining Ba — by. To avoid confusion it is worthwhile 
summarizing our choices. The quark mass m and Ca(^) are extracted simultaneously from 
the two-point AWL We then use these rh and c^(m) in all calculations of 6S. For improving 
the external current, Aj, in the three-point AWI we use ca{0). Lastly, the "slope-ratio" 
method, where we use the average of data with Cyi(m) and ca(0), gives cv{0) needed to 
improve the vector current. 

The estimate of b a — bp + bs/ 2 from fits to Eq. 1201 using the full set of masses (degenerate 
and non-degenerate) is very similar to that obtained using only the degenerate set. Including 
non-degenerate combinations we find that bm can also be extracted reliably. With bA — bp + 
bs/2 and bs in hand we can finally extract bp in two ways. The first is obtained by combining 
bp ~bs and bs and the second combines by, bA — by, bA — bp + bs/2 and 65. Both estimates 
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FIG. 17: Illustration of the signal for the ratio defined in Eq. (|18|) for the four data sets using 
two-point discretization. In all four cases the data have to be multiplied by the respective values 
of 2k3, the lattice normalization of the additional propagator in the numerator. 



use one combination of 6's from the three-point axial Ward identity. These are of similar 
quahty and dominate the errors. We find that these two estimates of bp, which provide a 
consistency check, differ at the level of the uncertainties present in all combinations of 6's 
extracted using three-point AWL For our final estimates of bp given in table EU we take the 
weighted average. 



XI. ct 

ct is extracted by solving, for each rhs, the Ward identity 

1 I ac ^.^<t-^4V-.](^^)(f,y4)Tif (0)) _ ^o E,-(^^r?^r(y^,y4)Tif (0)) 

and extrapolating these estimates to ms = as discussed in Ref. (1. The quality of the data 
for the ratios on the left and right hand sides of this equation is very good as illustrated 
in Figs. 1221 and 1221 We find that the two-point and three-point methods give consistent 
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FIG. 18: The 64NP two-point data for Zy/Z^Z^(m3) are obtained by extrapolating the ratio 
defined in Eq. ^Jto frii = m2 = using a quadratic fit. The intercept of the quadratic fit in m-a, 
gives Z^y/Z\Z\. 



estimates after the chiral extrapolations. We take the two-point value as our final estimate 
and the difference from the three-point result as a systematic uncertainty. 

The data, illustrated in Fig. |^ exhibit a behavior linear in ms. This can arise due to 
corrections of the form 0{aKam^). We had erroneously neglected this 0{m^a) dependence in 
ct in previous analyses. The slopes, —0.33(11), —0.21(10) and —0.17(3) at /3 = 6.0, 6.2 and 
6.4 respectively, are consistent with an aA behavior. The change in ct between a constant 
and a linear fit are significant at the la level, i.e. they change from 0.063(7) 0.085(12), 
0.051(7) ^ 0.063(10), 0.041(3) ^ 0.054(5) for the three (3 values respectively. Thus, our 
new estimates, based on linear fits, differ from those quoted in Ref. |^. 



XII. EQUATION-OF-MOTION OPERATORS 



We extract the coefficients, c^, of the equation-of-motion operators from the mi2 depen- 
dence of the three-point AWI 



7(13) 
'SO 



,(12) y{23) 



+ a^^±^mi2 + 0(a2). (22) 



'o 
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FIG. 19: is obtained from the product of ratios of correlators defined in Eqs. ()18() and ()13|) . The 
figure shows a quadric extrapolation in m3/2 for the 64NP data set with two-point discretization 
and c^(m). 



This can be rewritten as 

c'p + c'a = 2so - Xo (bso -bo- 26^) , (23) 

where Xo = Z^q/Z^ Zq and so is the slope, in the limit — 0, of the left hand side of 
Eq. (j22j) with respect to fhi for fixed m^. The results for dp + c'q are shown in Table EHl 
and for the three individual pieces so, Xoipso — bo)/'2-i and XohA in Table 1X1 

The quality of all the results is dominated by how well we can measure c'p. Unfortunately, 
the intermediate state in the relevant correlation functions is a scalar which has a poor signal. 
To obtain a fiat region with respect to the time slice of the operator insertion we fit the ratio 
on the l.h.s. of Eq. 1221 allowing m in 5S to be a free parameter. The resulting rh differ from 
those obtained using Eq. [7|by about 7%, 4%, and 2% at /3 = 6.0, 6.2 and 6.4 respectively. 

There is an additional systematic uncertainty of 0(a) ~ 0.1 in the determination of any 
slope from the chiral fits as discussed previously. This impacts the determination of all three 
terms so, Xoipso — bo)/ '2^ and XobA- 

Examples of fits to the left hand side of Eq. |221 are shown in Fig. |231 for the 64NP data 
set. The estimates, at leading order, should not depend on ma, however, the data show 
higher order effects. We, therefore, use a quadratic extrapolation in fh^ aX (3 = 6.4 and 
linear at /3 = 6.0 and 6.2 to get so- This changes the estimates from those presented in 
and 0. 



33 



-3.9 



-3 


95 




-4 


-4 


05 


-3 


98 




A 

— 4 


-4 


02 


-4 


04 


-4 


06 


-4 


04 


-4 


05 


-4 


06 


-4 


07 


-4 


08 


-4 


09 



^x60NPf 



^60NPb 5^ 

























— ^ 






LI 






■? 


V 
























^ < 


> < 


> <> 





























<> 







10 



20 



"1 1 1 r 



~i 1 1 r 



^ 62NP 



J I I L 



J I I L 



1 1 1 r 

J I I L 







10 



20 



30 



"1 1 1 r 



"1 1 1 r 



"1 1 1 r 



^ 64NP 



t I I I L 




J I I L 



' / \ ^ ^ 



J I I L 







10 



20 



30 



ay4 



FIG. 20: Comparison of the signal in the ratio of correlators on the l.h.s. of Eq. (|19|) used to extract 
Zp/Zg. In all four cases the data have to be multiplied by the respective values of 2k3, the lattice 
normalization of the additional propagator in the numerator. 



There is a very significant improvement in the signal for both the individual terms and the 
final c'p + c'q as [3 increases. Nevertheless, due to the uncertainties discussed above, all the 
c' could have additional systematic uncertainties similar to the errors quoted in Table IVllI, 
whose resolution is beyond the scope of this work. Thus, we consider our estimates as 
qualitative and warn the reader that the difference from the tree level value c'q = 1 should 
be used with caution. 



XIII. COMPARISON WITH RESULTS BY THE ALPHA COLLABORATION 

The ALPHA collaboration has used a very different method i.e., the Schrodinger Func- 
tional method, and their estimates have the largest differences from ours, so it is worthwhile 
comparing the two sets of values for Z^, Zy, Zp/Zg, ca, cy, and by- We expect the differ- 
ence to vanish as O(a^) for Z^, Zy and Zp/Z^, and as 0(a) for ca, cy and by. We find that, 
within combined errors, the estimates for by by the LANL, ALPHA and QCDSF collabora- 
tion are already consistent at all three (5 values as shown in Table HXl Similarly, estimates 
for Zp/Zg by the LANL and ALPHA collaborations agree. For the other four quantities, 
there is a statistically significant difference, and we have attempted to see whether the lattice 
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FIG. 21: Quadratic fits used to extract Z^Z^/Zp. The three fits correspond to (i) Eq. (|T9|) plotted 
versus afh = (ii) Eq. ()2U|) with fh defined using the mass dependent ca, and (iii) Eq. (|2U|) 

with m defined using the chirally extrapolated ca ■ The data are from the 64NP set with two-point 
discretization. Note that the intercepts from all three fits should agree up to errors of O(a^), but 
the slope of (i) is hp — bs whereas those of (ii) and (iii) give Ba — bp — bm- 

spacing dependence is consistent with theoretical expectations. To do this, we have fit the 
difference AX = X^anl — ^ alpha to an appropriate function of a, with the results: 

AZ^ = 0.004(1) - [261(16)a]2 X^/ndf = 0.03 (24) 

AZ^ = -0.002(12) + [222(190)a]2 x^/ndf = 0.5 (25) 

AcA = -188(39)a+ [769(74)a]2 xV^^^/ = 0.4 (26) 

Acy = -0.15(14) + 703(431)a x^/ndf = 0.01 (27) 

where a is in units of MeV~^ so that the coefficients are in units of MeV. Error estimates 
in AX were determined by adding the two independent statistical errors in quadratures. A 
number of comments are in order: 

• These fits are very sensitive to the errors assigned to AX and should only be used to 
draw qualitative conclusions. 
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FIG. 22: The signal in the ratio of correlators used to extract ct and defined on the left hand side 
of Eq. (|21() . The points are from 64NP two-point data with quark propagators. 

• As noted earlier, our condition [rriij independent of t for t >2) for fixing ca, and the 
variation in the physical size of our sources with a may lead to a more complicated 
dependence on a than simply the leading order expectation. Similarly, we have chosen 
different forms for the chiral extrapolation at the various /5's. These issues have been 
ignored here given the small number of values of (3. 

• For Zy and we expect a vanishing intercept and a difference proportional to a^. 
This expectation is borne out reasonably well. The non-zero value for the intercept in 
AZy could be a manifestation of higher order terms that are ignored in our fit. The 
size of the term is consistent with being ~ {qAqcd)^- 

• Fits to Aca without a quadratic term have large x^- The linear plus quadratic fit given 
in Eq. 123 does slightly better than constant plus quadratic. The fits are dominated by 
the difference at /5 = 6.0 where our estimate agrees with that given in Ref. |16i] . 

• Estimates of cy by the ALPHA collaboration are systematically much larger. The fit 
in Eq. |2ZI has large coefficients, however the errors are equally large. The calculation 
of cy warrants further study since the differences are large. 



36 



4.4 
4.3 



xeONPf 

oeoNPb <| >14-^^^ 



..<> 



<> 



<> 







10 



20 



4.4 F 
4.38 I 
4.36 
4.34 
4.32 



"1 1 1 r 



62NP 



"1 1 1 r 



T 1 1 r 



J I I L 



J I I L 



J I I L 







10 



20 



30 



4.36 



4.28 



"1 1 1 r 



"1 1 1 r 



4.34 

4.32 ^ 64NP 
4.3 



>( 



X 



T 1 1 r 



t I I I L 



J I I L 



J I I L 







10 



20 



30 



FIG. 23: The signal in the ratio of correlators defined on the right hand side of Eq. H21|) . This 
64NP data are used to extract ct- In all four cases the data have to be multiplied by the respective 
values of 2k3, the lattice normalization of the additional propagator in the numerator. 



XIV. COMPARISON WITH PERTURBATION THEORY 



The data at three values of the coupling allow us to also fit the difference between the 
non-perturbative and tadpole improved one-loop estimates as a function of a and a^, i.e., 
including both the leading order discretization and perturbative corrections. The results of 
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FIG. 24: ct is extracted using a quadratic fit to all seven quark masses (ki — ^7) for the 64NP 
data set with two-point discretization and Cyi(m). 



62iVP 


~\~ Cp 


so 


Xo{b5o-bo)/2 


XohA 


C'y + Cp 


-0.22(05) 


-0.06(2) 


1.50(5) 


+ ^'p 


-0.15(07) 


+0.05(5) 


1.42(5) 


Cp -\- Cp 


-0.17(10) 


-0.11(6) 


1.56(6) 


c'g + Cp 


-0.15(03) 


-0.05(1) 


1.29(5) 


~\~ Cp 


-0.21(07) 


+0.02(5) 


1.45(5) 


QANP 


c!q + c'p 


so 


Xo{h5o-ho)l2 


XohA 


c'y + dp 


-0.18(5) 


-0.08(3) 


1.36(6) 


+ Cp 


-0.08(6) 


+0.08(4) 


1.27(6) 


Cp -{- Cp 


-0.84(36) 


-0.41(32) 


1.49(8) 


C5 + Cp 


-0.06(4) 


-0.06(3) 


1.18(6) 


~\~ Cp 


-0.17(7) 


-0.05(5) 


1.31(6) 



TABLE X: The three contributions to the coefficient of the equation of motion operators c^ + Cp for 
the 62NP and 64NP data sets using the two-point derivative data and CA{fh) in the calculation 
of 5S and ca(0), cy(0), cr(0) in the discretization of the operators. 
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FIG. 25: Quadratic fits to the l.h.s. of Eq. versus am. The slopes so contribute to the 
coefficient of the equation of motion operators through Eq. (|23j) . The data are for 64NP with 
two-point discretization, CAifh) and ms chosen to be n^,- 

these fits are 



AZ^ = -{192af - {l.Sasf xV^df = 0.9 (28) 

AZ^ = -{159af - {1.2asf X^/ndf = 0.5 (29) 

A(Z°/Z°) = -(439a)2 - {Lda^f xVndf = 3.5 (30) 

Acy = -(138a) - (1.5a,)^ xVridf = 0.4 (31) 

Aca = (30a) - (lAasY X^/ndf = 0.01 (32) 

Act = (130a) + (OAasf X^/ndf = 0.1 (33) 

Abv = (1197a) - (S.Sasf X^/ndf = 3.2 (34) 

Abv = (630a) - (l.Sas)^ xVndf = 6.3 (35) 

AIa = -(770a) - {3.5asf xVridf = 0.8 (36) 

Abp = -(857a) - (1.9a,)2 x^/ndf = 2.3 (37) 

Abs = (507a) - (2.3a,)2 X^/ndf = 8.9 (38) 
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where a is expressed in MeV~^, AX = Xlanl — ^i-ioop, and as = g'^/ (47r'UQ) is the tadpole 
improved couphng with values 0.1340, 0.1255 and 0.1183 at the three (3. The tadpole factor 
uq is chosen to be the fourth root of the expectation value of the plaquette. 

One conclusion from these fits is that one-loop tadpole improved perturbation theory 
estimates of the Z's and c's underestimate the corrections. The deviations can, however, be 
explained by coefficients of reasonable size, i.e. the coefficient of 0(a) is ~ Aqcd and the 
perturbative corrections are (1 — 2)a^. The case of Zp/Zg is marginal, and we point to non- 
perturbative calculations using external quark and gluon states (the RI/MOM method) that 
show that the majority of the difference comes from 1-loop perturbation theory significantly 
underestimating (1 — Zp) [Ts^ . 

The most striking differences from perturbation theory are for the 6's. We stress, however, 
that the fits are very poor as evident from the /ndf . There are two useful statements we 
can make. In the case of hy (and similarly 6y), the agreement between our results and 
those by the ALPHA, QCDSF, and SPQcdR collaborations B H Q, suggests that 1- 
loop perturbation theory underestimates the correction. Second, at /5 = 6.4 6^, &p and 
hs are in good agreement with perturbation theory. 

XV. CONCLUSION 

We have presented new results for renormalization and improvement constants of bilinear 
operators at /? = 6.4. Combining these with our previous estimates at /3 = 6.0, and 6.2, and 
with the results from the ALPHA collaboration we are able to quantify residual discretization 
errors. Overall, we find that the efficacy of the method improves very noticeably with the 
coupling p. Using data at /3 = 6.4 we are able to resolve higher order mass dependent 
corrections in the chiral extrapolation for all the renormalization and improvement constants 
presented here. Our final results are summarized in Table IVll 

Determination of ca is central to 0(a) improved calculations. By comparing results 
from three different discretization schemes we improve the reliability of our error estimate. 
We also show that reliable estimates from correlators at finite momenta can be extracted 
and find that these give consistent results with those from zero-momentum correlators once 
additional O(p^a^) errors are taken into account. 

We find that both ca and cy are small, and the most significant differences from estimates 
by the ALHPA collaboration are at the strongest coupling (3 = 6.0. 

We also compare our non-perturbative estimates with one-loop tadpole improved pertur- 
bation theory. Overall, we find estimates based on 1-loop tadpole improved perturbation 
theory underestimate the corrections in the Z's and c's. The differences can, however, be 
explained by terms of 0{AqcdO') and (1 — 2)a^. 

The most significant differences are in by and by which are hard to explain by a com- 
bination of 0(a) and errors with coefficients of reasonable size. All the other &'s show 
agreement with perturbative estimates by /3 = 6.4. 
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